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COMPUTER SOLUTIONS OF THE GRAVITATIONAL N-BODY PROBLEM 


By Frank Hohl 

NASA Langley Research Center 

Thermalization effects in a one-dimensional model for a self -gravitating 
system have been investigated for systems containing low numbers of mass sheets 

(up to 4o) . It is found that the fluctuation of the kinetic energy is inversely 

proportional to the square root of the number of particles in the system. 

The stability of stationary states for one -dimensional self -gravitating 
systems is investigated by solving the N-body problem on an electronic com- 
puter. The number of "stars" in such systems is large (several 1000 ) so fiat 

the systems can be treated as collisionless . A variational method is used to 

show that stationary distribution functions that are always decreasing in going 
outward from the center of the system are stable. For other stationary distri- 
butions the system does not represent a minimum energy configuration and the 
system may be unstable . Computer solutions illustrating the resulting insta- 
bilities are presented. 

Numerical experiments with a simple two-dimensional rod model show that 
the spiral structure and other filamentary structure of galaxies may result 
from purely gravitational effects. 

THERMALIZATION EFFECTS IN A ONE-DIMENSIONAL STELLAR SYSTEM 


In a one -dimensional system crossings will always take place so that one 
can expect the system to approach a Maxwellian distribution. An exact double- 
precision computer program was used to investigate thermalization effects for 
low numbers of "stars." By an "exact" program we refer to a program where the 
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shortest crossing time for a mass sheet is determined and the system is then 
advanced by that time and the process is repeated. This, of course, differs 
from the model used for the results obtained for systems with large numbers of 
mass sheets where the system is always advanced for a constant time interval At, 
irrespective of the number of crossings during this time interval. The "exact" 
program is very accurate. For example, after time reversing the numerical inte- 
gration of a typical system investigated it was found that the initial condi- 
tions were reproduced accurately to 12 digits. 

The relaxation or characteristic time for a system of stars is given by 

T 0 = ( W )‘ l/2 ( 1 ) 

wheri - Ik the gravitational constant and p is the mass density. The 
characteristic length of interest is the Jeans or Debye length D which is 
defined by 

D = V T r c (2) 

where is the velocity dispersion of the system of stars. 

Lecar (ref. l) has investigated the exact one-dimensional motion of low 
numbers of "stars" and finds that to order hDt c there exists no thermaliza- 
tion. In the present study of thermalization effects a low number of mass 
sheets was chosen so that the system can be followed for a long period of time. 
In order to obtain meaningful results time averages of the systems are inves- 
tigated. The constant time interval used in the averaging process is taken to 
be smaller than the average crossing time for sheets in the system. The time 
averages are taken over times of the order n%)^T c where nD ~ N, the number 
of mass sheets in the system. It should be noted that the investigation of 



thermalization effects is still in progress and only some of the initial 
numerical results are presented here. 

For all of the systems investigated it was found that the time-averaged 
potential energy (P) and kinetic energy (T) satisfied the virial theorem. That 
is, <T>/<P> = 0.5 with an accuracy of at least three digits. The < > signify 
time-averaged quantities. As discussed by Ford (ref. 2) the Poincare recurrence 
time and the general behavior of the system is dependent on the initial con- 
ditions. For example, we can choose special initial conditions such that the 
time behavior of the system is very nearly periodic. However, for the results 
presented here the initial conditions were arbitrarily chosen. 

Figure 1 shows the time-averaged velocity distribution and density for a 
5-particle system. The solid line for the velocity distribution eorr'^’-onds 
to a Maxwellian distribution and is obtained by integrating A exp(-KU) over 
x. That is, 

f(v) = \\j^- | (5) 

The solid line for the density is given by 

n(x) = JtGm% 2 K sech^(jtGm%Kx) (4) 

The value of k is obtained from 


5N _ K 
2<T> + 2<E> 2<T> 


(5) 


The circles in figure 1 represent the time-averaged numerical results. The 
"experimental" velocity distribution is near the Maxwellian distribution. 
However, the variation of the density indicates that there is some near- 
periodic behavior of the system. Such near-periodic behavior is likely 
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to occur for very low numbers of "stars." Figure 2 shows the results for a 
four-sheet system. Again, the variation of the time-averaged velocity dis- 
tribution and density indicates that there is near-oscillatory behavior of the 
system. This oscillatory behavior has also been observed by plotting the 
kinetic energy of the system as a function of time. The fluctuations of the 
kinetic energy of the systems investigated are very violent and show no 
decrease in time. The fact that these fluctuations show no decrease in time 
indicates that a Maxwellian distribution can be reached while the fluctuations 
in the kinetic energy remain extremely violent. For example, figure 3 shows 
the time variation of the kinetic energy for a 6-particle system. The corre- 
sponding time-averaged velocity distribution and density are shown in figure 4. 
We see that the "experimental" points are near the Maxwellian distribution for 
both curves. Figure 5 shows that for a 10-particle system the "experimental" 
velocity distribution and density points reproduce the theoretical Maxwellian 
curves. The results presented so far were obtained by time averaging the sys- 
tem over times of the order of N^t c and the results indicate that the time- 
averaged distribution is very nearly Maxwellian. However, to show relaxation 
to a Maxwellian distribution one should prepare a system which is initially non- 
Maxwellian and then observe the relaxation to a Maxwellian distribution. This 
was attempted and the results for a 40-particle system are shown in figure 6. 

The solid line corresponds to a Maxwellian velocity distribution. The circles 
represent the initial velocity distribution of the system obtained by time 
averaging over 5 Nt c . After the system was advanced in time for 2N 2 charac- 
teristic periods it was found that the experimental velocity distribution was 
still identical to the initial distribution. Thus one must conclude that no 
relaxation to a Maxwellian distribution occurs to order W^t c . 
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Another quantity which gives information on the behavior of the system is 
the correlation function for the kinetic energy 

C(t) = <(T(t) - <C>)(T(t + t) - <!>)> (6) 


This quantity is plotted in figure 7 for a 6-particle system and the result 
shows that oscillations in the kinetic energy have a long-time memory and are 
only slowly damped. The energy correlation function for a 40-particle system 
is shown in figure 8. The results indicate that for small t there exists 
Landau damping hut some smaller correlation persists for a long time. If the 
time average is taken over a longer time, then the oscillations in C(t) are 
reduced for large r and remain the same for small t. In figure 9 "the rela- 
tive fluctuation of the kinetic energy 


d 


c(o) 

<T> 


(7) 


is plotted as a function of the number of particles N. It can be seen that 
the points fall near the curve 1 .fos which indicates that the oscillations in 
the kinetic energy of the systems do represent thermal fluctuations. It is 
therefore difficult to explain the long-time correlation shown in figures 7 
and 8. 


MINIMUM ENERGY PROPERTY 

It will now be shown that the minimum energy property which has been 
obtained by Hohl and Feix (ref. 3) for a special distribution can be extended 
to arbitrary distribution functions . The water-bag model illustrated in fig- 
ure 10 is used in the analysis. The contours v + ^)(x,t) and v_^^(x,t) 


describe surfaces of constant distribution function f = f^. According to the 
Liouville theorem the phase space bounded by the contours is incompressible. 
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In the limit of a very large number of contours the water-bag model can be used 
to construct arbitrary distribution functions. 

To simplify the equations we assume symmetric contours v*^) = v_(k) = v^) 
It is easily shown (ref. 4) that the equations which stationary contours v^) 
must satisfy are 

, X (k) 

v (k) dv E o (8) 

dx 


E is given by 


E(x) 



(9) 


where G is the gravitational constant, N is the number of mass sheets, each 
of mass m per -unit area, in the system. 

A e is defined by the distribution function 


f = £ aJs^v “ v - (k) ) - 5 ^v - v + ( k ^M 
k ~ . 

with 6_ 1 (z) = f Z The variable 9^) used in equation ( 9 ) is 


0 




v (k) U)d£ 


(10 ) 


(k) 

where x_ v ' ' is the end point of the. contour k. The- total energy of the sys- 
s 

tern is 
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Extremizing the integral for W requires that g satisfy the Euler-Lagrange 
equation 


or 







(k) dv 


00 


dx 


- E = 0 


( 12 ) 


(13) 


which are the equations for the stationary contours . 

If equations (13) are to represent a minimum energy configuration the 
Legendre's criterion of the second variation of g must he satisfied. That 
is, the quadratic form whose coefficient matrix has the elements 




d 2 g 


dv 




(J) 


(lit) 


must not he negative. Since only the diagonal elements 


a kk " 


00 


(15) 


are nonzero, the Legendre condition requires that 

A k 4o ( 16 ) 

for all k. Equation (l 6) is equivalent to stating that the distribution func- 
tion must always decrease in going outward from the center of the system where 
f = fj_ must he the largest. If equation (l6) is satisfied, the system is a 
minimum energy configuration and is always stable. However, if ^ 0 is not 
satisfied for all k, the system is not a minimum energy configuration. The 

• / V } 

system may then he unstable since the contours v v ' can now he deformed while 
keeping the total system energy constant. Numerical experiments with a one- 
dimensional model have been performed for two-contour systems to illustrate the 


7 



interchange instability which destroys the stationary state. For the two sys- 
tems investigated A-^ = -Ag . In the first case shown in figure 11 the ratio of 
the minimum to maximum star energy is 0.4. The figure shows that the stationary 
contours of the 2000-star system are quickly distorted while the heavier outer 
water bag tries to displace the inner hag. Figure 12 shows the results for a 
2000-star system with a ratio of minimum to maximum star energy equal to 0.25- 
The growth of the instability is now much slower because the central water bag 
or hole is much smaller. 

Another consequence of the minimum energy property is that any stable 
stationary state for which equation (l6) is satisfied can never be reached. 

TWO-DIMENSIONAL COMPUTER EXPERIMENTS FOR STELLAR SYSTEMS 

The motion of mass rods that are of infinite extent in the z-direction 
has been computed for 500-rod systems. The force acting on a particular mass 
rod is obtained by summing directly over the l/r force from each mass rod. 

This is a time-consuming process and the application of fast methods of solving 
the Poisson equation would speed up the calculations. A fast method of solving 
the Poisson equation has been used by Hockney (ref. 5) who followed the motion 
of 2000 mass rods . In comparing the results obtained by the two different 
methods it was found that they are nearly identical. 

The system of mass rods is advanced in time in the following manner. 

First, the force acting on all particles is computed by summing the l/r force 
for all particles. Second, the system is advanced for a small time step At 
and the process is repeated. The results of the calculations are displayed in 
x-y coordinate space. During the calculations the total energy and angular 
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momentum are computed to check on the accuracy of the computations . The nor- 
malizations 4jtG = 1 and m = 1 have been used for all the calculations. 

Figure 13 shows the time development of a system of 400 mass rods which 
has an initially rectangular distribution of uniform density in x-y space. 

The system has an initial thermal energy equal to l/5 of the initial potential 
energy plus an initial solid-body rotation equal to nearly twice that required 
to oppose the gravitational force toward the center of the system. It can be 
seen from figure 15 that the system quickly develops into a barred spiral. 
However, at a later time the spiral structure has almost completely disappeared 
and the system approaches a configuration similar to an elliptical galaxy. The 
time has been normalized to ca r “"*', the inverse of the frequency of the initial 
rotation . 

The remaining two-dimensional calculations were performed for 500 -particle 
systems which have an initially uniform circular distribution in x-y space 
and zero thermal velocity. The evolution of such systems is then studied for 
various values of inibial solid-body rotation. The initial positions are 
obtained by using a random number generator which gives a nearly -uniform dis- 
tribution over a circular region of the x-y plane. 

We first present the results for the case where the frequency of rotation, 
coy, equals cog, the frequency required such that the centrifugal force balances 
the gravitational force. Thus, 

% = cc g = j/hitGp (17) 

where p is the mass density of the rods per unit length. The resulting 
evolution of the system is shown in figure l4. The time has been normalized to 
ODg"-*-. Figure l4 shows that the system is relatively stable. At g = 6.32a^“- L 
there appear four irregular spiral arms. However, at a later time the spiral 
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arms almost completely disappear and the system takes an appearance reminiscent 
of an elliptical galaxy. 

The results for the case of zero initial rotation are presented next. Fig 
ure 15 shows that after an initial implosion the system expands and presents 
some highly irregular filamentary structure. After a second implosion the tem- 
perature of the system increases due to the randomness of the initial positions 
The pressure due to the temperature then tends to reduce the oscillations and 
the system again takes a form similar to an elliptical galaxy. 

For <%. = g- tDg the system again contracts initially and then expands . 

The results are shown in figure l6. An irregular structure appears initially 
which tends to disappear at a later time. Also at time t = 4.2a)g“^ the sys- 
tem is clustered into two aggregates which combine again at a later time. 

In figure 17 the results for the case cd*. = are shown. The system 

pulsates and shows some irregular structure. The general behavior is very 
similar to that of the previous case for o^. = The results for the two- 

dimensional stellar system are of a preliminary nature and additional work is 
required to investigate the evolution for systems with larger numbers of "stars 
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Kinetic energy 



Figure 5*- Time variation of the kinetic energy for a 6-particle system. 


















(b) Evolution up to t = J.i 
Figure 13.- Concluded. 



(a) Evolution up to t = 6.52o^“. 


Figure l4.- Evolution of a cylindrical stellar system with o^. = ccw 












(a) Evolution up to t = 1.86co g “ . 


Figure 15 •- Evolution of a cylindrical stellar system -with coj. 











(a) Evolution up to t = 3 . l 6 a 0 g -J - . 

Figure l 6 .- Evolution of a cylindrical stellar system with ooj. 








(a) Evolution up to t = 6.32cBg-- L . 

Figure IT.- Evolution of a cylindrical stellar system with 








(t>) Evolution up to t = 


Figure 1J,- Conclud 
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